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ABSTRACT 

Context. Detection of solar gravity modes remains a major challenge to our understanding of the inner parts of the Sun. Their 
. frequencies would enable the derivation of constraints on the core physical properties, while their amplitudes can put severe 

■ constraints on the properties of the inner convective region. 

, Aims. Our purpose is to determine accurate theoretical amplitudes of solar g modes and estimate the SOHO observation 

■ duration for an unambiguous detection of individual modes. We also explain differences in theoretical amplitudes derived from 
' previous works. 

, Methods. We investigate the stochastic excitation of modes by turbulent convection, as well as their damping. Input from a 
' 3D global simulation of the solar convective zone is used for the kinetic turbulent energy spectrum. Damping is computed 
using a parametric description of the nonlocal, time-dependent, convection-pulsation interaction. We then provide a theoretical 
I estimation of the intrinsic, as well as apparent, surface velocity. 

' Results. Asymptotic g-mode velocity amplitudes are found to be orders of magnitude higher than previous works. Using 
a 3D numerical simulation from the ASH code, we attribute this to the temporal-correlation between the modes and the 

■ turbulent eddies, which is found to follow a Lorentzian law rather than a Gaussian one, as previously used. We also find 
' that damping rates of asymptotic gravity modes are dominated by radiative losses, with a typical life time of 3 x 10^ years 
, for the i — 1 mode ai f = 60/iHz. The maximum velocity in the considered frequency range (10-100 ^Hz) is obtained for 

■ the i = 1 mode at = 60 ^Hz and for the £ = 2 at = 100 ^Hz. Due to uncertainties in the modeling, amplitudes at 
maximum i.e. for £ = 1 at 60 /iRz can range from 3 to 6 mm s~^. The upper limit is too high, as g modes would have been 

. easily detected with SOHO, the GOLF instrument, and this sets an upper constraint mainly on the convective velocity in the Sun. 



Key words, convection - turbulence - Sun: oscillations 



1. Introduction 

The pioneer works of Ulrich (1970) and Leibacher & Stein 
(1971) led to the identification of the solar five-minutes 
oscillations as global acoustic standing waves (p modes). 
Since then, successful works have determined the Sun in- 
ternal structure from the knowledge of its oscillation fre- 
quencies (e.g., Christensen-Dalsgaard 2004). However, p 
modes are not well-suited to probing the deepest inner 
part of the Sun. On the other hand, g modes are mainly 
trapped in the radiative region and are thus able to pro- 
vide information on the properties of the central part of 
the Sun (r < 0.3 i?©) (e.g., Turck-Chieze et al. 2001; 
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Christensen-Dalsgaard 2006). As g modes are evanescent 
in the convective region, their amplitudes are expected to 
be very low at the photosphere and above, where obser- 
vations arc made, making their detection is thus quite a 
challenge for more than 30 years. 

The first claims of detection of solar gravity modes be- 
gan with the work of Severnyi et al. (1976) and Brookes 
et al. (1976). Even after more than ten years of observa- 
tions from SOHO, there is still no consensus about detec- 
tion of solar g modes. Most of the observational efforts 
have been focused on low-order g modes motivated by 
a low the granulation noise (Appourchaux et al. 2006; 
Elsworth et al. 2006) and by previous theoretical esti- 
mates of g-mode amplitudes (e.g., Turck-Chiczc et al. 
2004; Kumar ct al. 1996). Recently, Garcia et al. (2007) 
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have investigated the low-frcqucncy domain, with the hope 
of detecting high radial-order g modes. The method looked 
for regularities in the power spectrum, and the authors 
claim to detect a periodicity in accordance with what is ex- 
pected from simulated power spectra. The work of Garcia 
et al. (2007) present the advantage of exploring a different 
frequency domain {i/ € [25; 140]/iHz) more favorable to a 
reliable theoretical estimation of the g-mode amplitudes, 
as we will explain later on. 

Amplitudes of g modes, as p modes, are believed to 
result from a balance between driving and damping pro- 
cesses in the solar convection zone. Two major processes 
have been identified as stochastically driving the reso- 
nant modes in the stellar cavity. The first is related to 
the Reynolds stress tensor, the second is caused by the 
advection of turbulent fluctuations of entropy by turbu- 
lent motions. Theoretical estimations based on stochas- 
tic excitation have been previously obtained by Gough 
(1985) and Kumar et al. (1996). Gough (1985) made an 
order of magnitude estimate based on an assumption of 
equipartition of energy as proposed by Goldreich & Keeley 
(1977b). He found a maximum velocity around 0.5 mms~^ 
for an i = 1 mode at w 100 /xHz. Kumar et al. (1996) 
used a different approach based on the Goldreich et al. 
(1994) modeling of stochastic excitation by turbulent con- 
vection, as well as an estimating of the damping rates 
(Goldreich & Kumar 1991) that led to a surface velocity 
near 0.01 mm s~^ for the £ = 1 mode at « 100 /iHz. 
The results differ from each other by orders of magni- 
tude, as pointed out by Christenscn-Dalsgaard (2002b). 
Such differences remain to be understood. One purpose of 
the present work is to carry out a comprehensive study 
of both the excitation and damping rates of asymptotic g 
modes. Our second goal is to provide theoretical oscillation 
mode velocities, , as reliably as possible. Note, however, 
that penetrative convection is another possible excitation 
mechanism (Andersen 1996; Dintrans et al. 2005), but it 
is beyond the scope of this paper. 

Damping rates arc computed using the Grigahccne 
et al. (2005) formalism that is based on a non-local time- 
dependent treatment of convection. We will show that, 
contrary to p modes and high frequency g modes, asymp- 
totic g-mode (i.e low frequency) damping rates are insen- 
sitive to the treatment of convection. This then removes 
most of the uncertainties in the estimated theoretical os- 
cillation mode velocities. Consequently, we restrict our in- 
vestigation to low-frequency gravity modes. Stochastic ex- 
citation is modeled in the same way as in Belkacem et al. 
(2008), which is a generalization to non-radial modes of 
the formalism developed by Samadi & Goupil (2001) and 
Samadi et al. (2003b, a), for radial modes. As in the case 
of p-modes, the excitation formalism requires knowing the 
turbulent properties of the convection zone, but unlike p 
modes, the excitation of gravity modes is not concentrated 
towards the uppermost surface layers. One must then have 
some notion about the turbulent properties across the 
whole convection zone. Those properties will be inferred 



from a 3-D numerical simulation provided by the ASH 
code (Miesch et al. 2008). 

The paper is organized as follows. Section 2 briefly re- 
calls our model for the excitation by turbulent convection 
and describes the input from a 3D numerical simulation. 
Section 3 explains how the damping rates arc computed. 
Section 4 gives our theoretical results on the surface veloc- 
ities of asymptotic g modes and compares them with those 
from previous works. Section 5 provides the apparent sur- 
face velocities, which take disk integrated effects and line 
formation height into account. These quantities can be 
directly compared with observations. We then discuss our 
ability to detect these modes using data from the GOLF 
instrument onboard SOHO as a function of the observing 
duration. The discussion is based on estimations of detec- 
tion threshold and numerical simulations of power spectra. 
In Sect. 6., uncertainties on the estimated theoretical and 
apparent velocities, due to the main uncertainties in our 
modeling, are discussed. Finally, conclusions are provided 
in Sect. 7. 



2. Excitation by turbulent convection 

The formalism we used to compute excitation rates of non- 
radial modes was developed by Belkacem et al. (2008) who 
extended the work of Samadi & Goupil (2001) developed 
for radial modes to non-radial modes. It takes the two 
sources into account that drive the resonant modes of the 
stellar cavity. The first is related to the Reynolds stress 
tensor and the second one is caused by the advection of 
the turbulent fluctuations of entropy by the turbulent mo- 
tions (the "entropy source term" ). Unlike for p modes, the 
entropy source term is negligible for g modes. We numer- 
ically verified that it is two to four orders of magnitude 
lower than the Reynolds stress contribution depending on 
frequency. This is explained by the entropy contribution 
being sensitive to second-order derivatives of the displace- 
ment eigentunctions in the superadiabatic region where 
entropy fiuctuations arc localized. As the gravity modes 
are evanescent in the convection zone, the second deriva- 
tives of displacement eigentunctions arc negligible and so 
is the entropy contribution. 

The excitation rate, P, then arises from the Reynolds 
stresses and can be written as (see Eq. (21) of Belkacem 
et al. 2008) 

3 rM r+oo 

P=— dmpnR{r) / dkSk (1) 
Jo Jo 

Sk^^ j du (k) Xkiio + u;q) Xkiio) (2) 

where m is the local mass, po the mean density, luq the 
mode angular frequency, / the mode inertia, Sk the source 
function, E{k) the spatial kinetic energy spectrum, Xk the 
eddy-time correlation function, and k the wavenumber. 
The term R{r) depends on the eigenfunction, its expres- 
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sion is given in Eq. (23) of Belkacem ct al. (2008), i.e 
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and are the radial and horizontal components of the 

fluid displacement eigenfunction (^), and £,m represent 
the degree and azimuthal number of the associated spher- 
ical harmonics. 
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Fig. 1. Luminosity flux contributions versus radius, averaged 
over horizontal surfaces and in time. The solid line corresponds 
to the enthalpy luminosity (Le), the short dashed- line to the 
radiative luminosity (Lr), the dotted-dashed line to the ki- 
netic energy luminosity ( Luin), the long-dashed line to the 
total luminosity (Ltot), and the dot-dashed line correspond to 
the unresolved eddy luminosity (Led) (Brun et al. 2004). We 
particularly emphasize the negative kinetic-energy flux that 
results in a larger convective flux (see text for details). 



2.1. Numerical computation of theoretical excitation 
rates 

In the following, we compute the excitation rates of 
g modes for a solar model. The rate {P) at which en- 
ergy is injected into a mode per unit time is calculated 
according to Eq. (1). Eigenfrequencies and eigenfunctions 
are computed using the adiabatic pulsation code OSC 
(Boury et al. 1975). The solar structure model used for 
these computations is obtained with the stellar evolu- 
tion code CESAM (Morel 1997) for the interior and a 
Kurucz (1993) model for the atmosphere. The interior- 
atmosphere matching point is chosen at logr = 0.1 
(above the convective envelope). The pulsation compu- 
tations use the full model (interior-l- atmosphere) . In the 
interior model, we used the OPAL opacities (Iglesias & 
Rogers 1996) extended to low temperatures with the opac- 
ities of Alexander & Ferguson (1994) and the CEFF equa- 
tion of state (Christensen-Dalsgaard & Dappen 1992). 
Convection is included according to a Bohm-Vitense 
mixing-length (MLT) formalism (see Samadi et al. 2006, 
for details), from which the convective velocity is com- 
puted. Turbulent pressure is not included (but see discus- 
sion in Sect. 6). 

Apart from the eigenfunctions and the density strati- 
fication, Eq. (1) involves both the convective velocity and 
the turbulent kinetic energy spectrum. To get some in- 
sight into the turbulent properties of the inner part of 
the solar convection zone, we chose to use results from 
(ASH) 3D numerical simulations. Such a choice was mo- 
tivated by the uncertainties inherent in the treatment of 
turbulence by the MLT. The MLT indeed only gives us 
an estimation of the convective flux but is not able to as- 
sess the contributions of all scales involved in turbulent 



convection. Thus, in the following, the rms convective ve- 
locity is taken from the mixing-length theory, while both 
the spatial and temporal turbulent properties are inferred 
from the 3D simulation. Then, velocity from the numeri- 
cal simulation is not used in our calculation. This choice 
is motivated by the rigid boundary condition at the top 
of the simulation that results in an unrealistic decrease in 
the vertical velocity for r > 0.93 i?©. 

2.2. The 3D convection simulation 

One way of assessing the dynamical properties of the deep 
solar turbulent convection zone is to exploit a high reso- 
lution numerical simulation such as those performed with 
the anelastic spherical harmonic (ASH) code (Miesch et al. 
2008; Brun et al. 2004). The simulation of global scale tur- 
bulent convection used in the present work is discussed in 
detail in Miesch et al. (2008). The ASH code solves the 
hydrodynamic anelastic equations within a spherical shell 
extending from r = 0.71uptor = O.98i?0, yielding an 
overall radial density contrast of 132. Solar values were 
assumed for the rotation rate and the imposed luminos- 
ity. Figure 1 represents the energy flux balance (converted 
to luminosity and normalized to the solar luminosity) in 
the simulation. We clearly see how dominant, and overlu- 
minous, the convective (enthalpy) flux is in carrying the 
heat outward. This is mostly due to the strong density 
contrast and to the corresponding strong asymmetry be- 
tween up- and downflows yielding a large inward kinetic 
energy flux (see Miesch et al. 2008 for more details). We 
have seen above that, in order to compute the excitation 
rate of the waves, one needs some well-defined physical 
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Fig. 2. E{kh) computed as explained in Appendix. B, for three 
shell radii that sample the convection zone, as a function of the 
local horizontal wave number kh- 



quantities, such as the kinetic energy spectrum {Ek) and 
the eddy time function {xk)- It is straightforward to de- 
duce these quantities from the 3D simulation as explained 
in Appendix B. We then directly use Ek in Eq. (2) to 
compute the source function, whereas for Xk we perform 
a fit of the 3-D results with a simple analytical expres- 
sion. In the ASH code, the set of anelastic equations is 
projected onto spherical harmonics for the horizontal di- 
mensions. This implies that the kinetic energy spectrum 
is obtained as a function of the spherical degree / . The 
local wavenumber kh is obtained via the simple expression 
kh = \/K^ + "^ith r the shell radius. 

2.2.1. Kinetic energy spectrum and time-correlation 
function 

The kinetic energy spectrum of the total velocity {i.e. the 
horizontal and vertical components), E{kh), is plotted in 
the top panel of Fig. 2 as a function of the local horizon- 
tal wave number kh. The rms convective velocity (it) in- 
creases with r, thus explaining that the deeper the layers, 
the smaller E{k) since J dkE{k) ~ 1/2 w^. In terms of ex- 
citation rates, an important issue is the scale at which the 
spectrum peaks. As pointed out by Miesch et al. (2008), 
the scale at which the kinetic energy spectrum is maxi- 
mum is the scale between the downwflows. It is about 58 
Mm at the top of the simulation (r — 0.98 i?©) up to 300 
Mm at the bottom. This is quite different from what is 
found in the uppermost layers in 3D numerical simula- 
tions of the Sun (e.g.. Stein & Nordlund 1998), in which 
the maximum of E{k) is found on a scale around 1 Mm. 
Such a difference is explained by the density that strongly 
decreases in the upper layers. 

The time-correlation function (xk) also plays an im- 
portant role. Usually, a Gaussian time-correlation function 
is used (Goldreich et al. 1994; Chaplin et al. 2005). Samadi 
et al. (2003a) demonstrate that Xk is reproduced better by 
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Fig. 3. Top: Crosses represent Xfe('^) obtained from the 3D 
simulation at the wave number fco that corresponds to the 
maximum of E(k), and at the radius t/Rq = 0.89. Data are 
obtained with a time series of duration « 45.83 days with a 
sampling time of 4 • 10^ seconds. Analytical functions are nor- 
malized so that their integrals are equal to unity. Bottom: 
The same as the upper panel except that data are obtained 
with a time series of duration ~ 4.68 days with a sampling 
time of 800 seconds. The theoretical curves are normalized so 
that their integrals over frequency equal that of the simulated 
data. 



a Lorentzian function. They argue that the departure from 
a Gaussian function can be explained by the presence of 
plumes in the uppermost part of the convection zone. This 
result, obtained with 3-D numerical simulations, was then 
confirmed by confronting solar-p modes excitation rates, 
computed with Gaussian and Lorentzian functions, with 
the observational data. It turns out that the Lorentzian 
function greatly improves the agreement between models 
and observations. However, the time-correlation function 
is unknown at deeper layers. The eddy-time correlation 
function derived from the 3D numerical simulation pro- 
vided by the ASH code is therefore compared to Gaussian 
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and Lorcntzian functions that arc respectively defined as 
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where ujk is its linewidth, defined as 
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where A is a parameter as in Balmforth (1992). the velocity 
Ufc of the eddy with wavenumber k is related to the kinetic 
energy spectrum E{kh) by (Stein 1967) 



2k 



dkE{k) 



(11) 



Figure 3 presents the comparison between analyti- 
cal time-correlation functions, computed following the 
set of Eqs. (7)-(ll), and Xk computed from the 3D nu- 
merical simulation. The latter is calculated as described 
in Appendix. B. The Lorentzian function represents the 
eddy-time correlation function better than a Gaussian 
function in the frequency range we are interested in (i/ e 
[20/iHz;110AiHz]). 

The best fit is found using a sum of a Lorentzian func- 
tion with A = 3 and a Gaussian with A = 1/3 as shown in 
the top panel of Fig. 3. In the frequency range we are in- 
terested in, i.e. at frequencies corresponding to the gravity 
modes (bottom panel of Fig. 3) the fit reproduces the time- 
correlation given by the 3-D numerical simulation. We also 
clearly see that the eddy-time correlation function is very 
poorly represented by a Gaussian function, which only 
reproduces very low frequencies that do not significantly 
contribute to the excitation, then it fails and underesti- 
mates Xk by many order of magnitudes (see Sect. 4.2.2). 

The results presented in Fig. 3 are for the depth 
r Ki O.SRq, where excitation is dominant, and for an angu- 
lar degree corresponding to the maximum of the kinetic 
energy spectrum [i = 40), whose contribution is domi- 
nant in the excitation rates. Those results do not depend 
on the shell considered but instead on the wavenumber. 
For very high angular degree {£ > 300) we find that Xk 
becomes more and more Gaussian. Nevertheless, as shown 
by Fig. 2, those contributions are negligible compared to 
large-scale ones. 

The value of the parameter A is also of interest. 
Contrary to the upperlayers where A = 1 (Samadi et al. 
2003b), we find a higher value, A = 3, that accords with 
the result of Samadi et al. (2003b) who find that the 
deeper the layers, the higher this parameter. 




d5(\,Tj) at r/Ro = 0.74- 




Fig. 4. The source function is plotted versus the spherical an- 
gular degree (/), and the frequency for two radii: r — 0.957?© 
(top panel) and r — O.747?0 (bottom panel). Bright (red) and 
dark (blue) tones indicate the high and low intensity of the 
source function, respectively. The color table is logarithmic. 
The black line corresponds to an arbitrary contour line that is 
the same for both panels. 



2.2.2. The source function 

Figure 4 displays the source function {Sk, Eq. (1)) as a 
function of both the angular degree I involved in the sum- 
mation Eq. (B.l) and the mode frequency. The function 
Sk evaluated at two levels, r = O.95i?0 and r = 0.747?©, 
is shown in order to emphasize the dependence of Sk with 
the radius. Near the top of the convection zone, Sk is non- 
negligible at high frequencies (i/ > 50/iHz) and on small 
scales. From top to bottom, the intensity of the source 
function decreases such that at the bottom, significant in- 
tensities exist only on large scales (small I values) and low 
frequencies. This behavior corresponds to the evolution of 
convective elements, i.e turbulent eddies evolve on larger 
time and spatial scales with depth. Thus, we conclude that 
high-frequency g modes are mainly excited in the upper 
layers, whereas low ones are excited deeper, however, the 
net excitation rate, Eq. (1). is a balance between the eigen- 
function shape and the source function. 
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2.3. Excitation rates 



Anticipating the following (see Sect. 3), we stress that 
modes with high angular degree will be highly damped, 
making their amplitudes very small; hence, we restrict our 
investigation to low-^ degrees {i < 4). In Fig 5, we present 
the excitation rates for low-frequency gravity modes (i.e., 
£=1,2,3). By asymptotic modes we denote low- frequency 
modes (i^ < 100 /^Hz, i.e high-|n| modes) while high fre- 
quencies {v > 100 /uHz) correspond to low-|n| modes. At 
low frequencies (z/ < 100 /xHz), the excitation rate (P) de- 
creases with increasing v, it reaches a minimum and then 
at high frequency increases with the frequency. This can 
be explained by considering the two major contributions 
to the excitation rate P (Eq. (1) and Eq. (3)) which are 
the inertia / (in Eq. (1)) and mode compressibility (V • ^, 
appearing in i?(r), Eq. (3)). 

Mode inertia decreases with frequency as shown by 
Fig. 6 since the higher the frequency, the higher up the 
mode is confined in the upper layers. This then tends to 
decrease the efficiency of the excitation of low-frequency 
modes. On the other hand, mode compressibility (Fig. 6) 
increases with frequency and consequently competes and 
dominates the effect of mode inertia. Mode compressibility 
can be estimated as 



dr 



(12) 



The mode compressibility is minimum when both terms in 
Eq. (12) are of the same order. Following Belkacem et al. 
(2008), one has 



d^r /(£+l)a 



dr 



with 



(13) 



where a is the dimensionless frequency, ojq is the angular 
frequency of the mode, R the Sun radius, and M its mass. 
According to Eq. (13), mode compressibility is minimum 
for V K, 100 /^Hz depending on as shown by Fig. 6. In 
contrast, in the asymptotic regime {v < 100 /iHz), the 
modes are compressible and this compressibility increases 
with decreasing frequency. 

It is important to stress that for the asymptotic 
g modes, in the frequency range [20; 110] /xHz, the hor- 
izontal contributions in Eq. (3) are dominant. For low-£ 
g modes, the dominant contributions come, in Eq. (3), 
from the component of the mode divergence (see Eq. (12)). 
Then the ratio of the horizontal to the vertical contribu- 
tions to Eq. (1) is around a factor five, imposing the use 
of a non-radial formalism. 

3. Damping rates 

To compute theoretical (surface velocities) amplitudes of 
g modes, knowledge of the damping rates is required. 



3.1. Physical input 

Damping rates have been computed with the non- 
adiabatic pulsation code MAD (Dupret 2002). This code 
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Fig. 5. Rate (P) at which energy is supplied to the modes 
versus the frequency for modes with angular degree £ = 1,2, 
and 3. The computation is performed as detailed in Sect. 2.1, 
using a Lorentzian eddy-time correlation function. 



includes a time-dependent convection (TDC) treatment 
described in Grigahcene et al. (2005): it takes into ac- 
count the role played by the variations of the convcctive 
flux, the turbulent pressure, and the dissipation rate of 
turbulent kinetic energy. This TDC treatment is non-local, 
with three free parameters a, 5, and c corresponding to the 
non-locality of the convcctive flux, the turbulent pressure 
and the entropy gradient. We take here the values a = 10, 
6 = 3, and c = 3.5 obtained by fitting the convcctive fiux 
and turbulent pressure of 3D hydrodynamic simulations 
in the upper overshooting region of the Sun (Dupret et al. 
2006c). According to Grigahcene et al. (2005), we intro- 
duced a free complex parameter /3 in the perturbation of 
the energy closure equation. This parameter is introduced 
to prevent non-physical spatial oscillation of the cigen- 
functions. We use here the value /3 = —0.5i, which leads 
to a good agreement between the theoretical and observed 
damping rates and phase lags in the range of solar pres- 
sure modes (Dupret et al. 2006a). The sensitivity of the 
damping rates to f3 is discussed in Sect. 3.2.1, and we 
show in next sections that the values of those parameters 
have no influence on the results since we are interested in 
low-frequency g modes. 

We use the TDC treatment as described in Dupret 
et al. (2006b), in which the ID model reproduces exactly 
the mean convcctive flux, the turbulent pressure and the 
mean superadiabatic gradient obtained from a 3D hydro- 
dynamic simulation by Stein & Nordlund (1998), by in- 
troducing two fitting parameters, the mixing-length, and 
a closure parameter (see Dupret et al. 2006b, for details). 
We also stress that, for low-frequency g modes, particular 
attention is to be paid to the solution of the energy equa- 
tion near the center as explained in Appendix A for the 
£ = I modes since those dipolar modes present a peculiar 
behavior near the center that must be properly treated. 
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Fig. 6. Top: Absolute value of mode compressibility for 1= \ 
modes versus the frequency, computed for three different lay- 
ers in the convection zone. Bottom: Mode inertia versus fre- 
quency for modes with angular degree i= 1, 2, 3. 
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Fig. 7. Work integrals for i= \ and 1 = 2 modes at ~ 60/iHz 
(top panel) and v ~ 20/iHz (bottom panel), the surface values 
give the damping rates r\ in /xHz. 



3.2. Numerical results for a solar model 

3.2.1. Sensibility to the time-dependent treatment of 
convection 

To understand the contribution of each layer of the Sun in 
the damping of the g modes, we give the normalized work 
integral in Fig. 7 in such a way that the surface value is 
the damping rate 77 (in /xHz)^. Results obtained with our 
TDC treatment (solid lines) and with frozen convection 
(FC, dashed line) are compared for 4 different modes with 
V ~ 60/xHz (top panel) and v ~ 20/.tHz (bottom panel). 
We see that most of the damping occurs in the inner part 
of the radiative core. The work integrals obtained with 
TDC and FC treatments are not very different; hence, 
the uncertainties inherent in the treatment of the coher- 
ent interaction between convection and oscillations do not 
significantly affect the theoretical damping rates of solar 
asymptotic g modes. This means that the frozen convec- 
tion is adapted to low-frequency g modes. This can be 

^ Note that, regions where the work decreases outwards have 
a damping effect on the mode, or a driving effect when it in- 
creases outwards. 



explained by paying attention to the ratio Q = Wg/wc, 
where loq is the oscillation frequency and ujc the convec- 
tive frequency, defined to be ujc = 27rA/u,„;t where A is 
the mixing length and Umit the convective velocity. In the 
whole solar convective zone Q is higher than unity ex- 
cept near the surface (the superadiabatic region). However 
contributions of the surface layer remain small in compar- 
ison with the radiative ones for asymptotic g modes (see 
Fig. 7). 

One can thus draw some conclusions 

— for high-frequency g modes {v > 110 /xHz), the work 
integrals and thus the damping rates are sensitive to 
the parameter /3 that is introduced to model the con- 
vection/pulsation interactions because the role of the 
surface layers in the work integrals becomes impor- 
tant. As a result, the results on the damping rates are 
questionable for high frequencies since the value of (3 
is derived from the observed p modes and that there 
is no evidence it can be applied safely for g modes. 

— in contrast, for low-frequency g modes {v < 110 /iHz), 
we find that the work integrals and then the damping 
rates are insensitive to parameter (3. Also, we numeri- 
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Fig. 8. Contributions to the work by the radial radiative 
flux variation (solid line), the transverse radiative flux vari- 
ation (dotted line), and the time-dependent convection terms 
(dashed line), for the mode 1 = 1, gio (top panel) and £ — 1, 
g32 (bottom panel). Details are given in the text. 

cally checked that the damping rates are insensitive to 
the non-local parameters introduced in Sect. 3.1. 



3.2.2. Contributions to the work integral 

Figure 8 allows us to investigate the respective roles played 
by different terms in the damping of the mode. More pre- 
cisely we consider two modes {i = 1, gio, v « 60 /iHz 
and i = 1, g32 , ~ 20 ^Hz) in the frequency interval of 
interest here and give in Fig. 8 the modulus of: 

— the contribution to the work by the radial part of the 
radiative flux divergence variations (solid line) 



FRr 



5T 



dSL 



R 



R 



dx ' GAPa 



(14) 



where T is the temperature, Lj^ the radiative lumi- 
nosity, R the solar radius, M the solar mass, x the 
normalized radius, a the real part of the normalized 
frequency a ~ ujo/ [GM / R^Y^ , and x the normahzed 



10" 



10" 
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Fig. 9. Theoretical damping rates of g modes of degree I — 
1, 2, 3 as a function of the oscillation frequency in ^Hz. 

radius (see Appendix A). Note that S denotes the wave 
Lagrangian perturbations, 3? the real part, and * the 
complex conjugate. 
— the contribution to the work by the transversal part 
of the radiative flux divergence variations (dotted line): 



FRh 



-^(£+1)3? 



5T* 



5T 



r;dT/dx 



RL 



(15) 

— the contribution to the work by the time-dependent 
convection terms (dashed line): dWc (see Sect. 4 of 
Grigahcene et al. 2005). 

Integration of these terms over the normalized radius gives 
their global contribution to the work performed during one 
pulsation cycle. 

The time-dependent convection terms have a very low 
weight for both modes in the frequency range ly < llO yuHz. 
It confirms the conclusion of Sect. 3.2.1 that the damping 
rates of low-frequency g modes are not dominated by the 
perturbation of the convcctive flux, i.e. the interaction con- 
vection/oscillation (through the parameter (3). The higher 
the mode frequency, the higher the integrated convective 
contribution of the work (Wc), which becomes dominant 
for > 110 fiRz. 

While the transverse radiative flux term plays a signif- 
icant role near the center, the major contribution to the 
work comes from the radial component of the radiative 
flux variations. As a result, the radiative damping is the 
dominant contribution for low- frequency gravity modes. 

In Fig. 9, we give the theoretical damping rates rj of g- 
modes of degree £ = 1,2,3, as a function of the oscillation 
frequency in fiHz. We sec that for i/ < 110 /iHz, 77 is a de- 
creasing function of frequency. We find that the frequency 
dependence is cx cjq^. To understand this behavior, we 
express the integral expression of the damping rate (see 
Grigahcene et al. 2005, for details) as 



J— r im (^-^'t5s] (Fs - l)dm 



(16) 
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with 



dm lip and (Fg - 1) = 



fdlnT\ 



(17) 



where Sp, SS are the perturbations of the density and en- 
tropy, respectively, p, T are the density and temperature, 
I the cigcnfunction, and the star denotes the complex con- 
jugate. 

Keeping only the radial contribution of the radiative 
flux in the energy equation (Eq. (A.l) ) because it is the 
dominant contribution, and neglecting the production of 
nuclear energy (e = 0), one gets 



TSS 



i dSL 
luq dm 



(18) 



This approximation comes from the dominance of the ra- 
dial contribution of the radiative flux. In addition, in the 
diffusion approximation 



1 



dST 



ST 
~T 



6k 6p d^r 
K p dr 



5L _ I 

T ^ \ (dT/dr) 

(19) 

Because of the high wavenumber for low-frequency g 
modes, the term in dST/dr is very high in Eq. (14), domi- 
nates in Eq. (19), and is the main source of damping. This 
term appears as a second-order derivative in the work inte- 
gral, and introduces a factor fc^ (fc^ ~ y/J(T^rT)N/ (tuor) is 
the vertical local wavenumber). Thus, from Eqs. (19), (18), 
and (16) one obtains rj oc ujq^/I. By using an asymptotic 
expansion of the eigenfunctions (Christensen-Dalsgaard 
2002a), one gets / cx liJq'^, which permits rj cx ujq^ and 
explains the behavior of 77 in Fig. 9. The argument is the 
same for the variation of 77 with the angular degree at 
fixed frequency becomes it comes from the wave-number 
dependence kf.. 

Above 110 /iHz, the role of the radiative zone in the 
mode damping is smaller. There, the damping rates be- 
gin to increase with frequency simply because the kinetic 
energy of the modes decreases faster than the mechanical 
work. 



4. Surface velocities of g modes 

4.1. Theoretical (intrinsic) velocities 

We compute the mean-squared surface velocity (w^) for 
each mode as 

v'sih) = Q^^(v{r,t)-vir,t))dnj [h) (20) 

where h is the height in the stellar atmosphere, <> the 
time average. Using the expression Eq. (C.3) in appendix 
C, one then has 

vlih) = [vlih) + t{i + l)vl{h)\ . (21) 

The amplitude = (1/2) {\a{t)\'^) is given by (Eq. (C.6)): 



(22) 



where <> denotes the time average, / the mode inertia, 77 
the damping rate, and Vj,,hih) = luq ^r.h{h) with ^r(^) and 
^hih) respectively the radial and horizontal displacement 
eigenmode components. 

In this section, we consider the level of the photo- 
sphere h ^ R with R the radius at the photosphere. 
Figure 10 presents intrinsic values of the velocities. The 
behavior of the surface velocities as a function of the 
angular degree {£) is mainly due to the damping rates, 
which rapidly increase with £; hence, at fixed frequency, 
the higher the angular degree, the lower the surface veloc- 
ities. As a consequence, amplitudes are very low for £ > 3 . 
At fixed £, Vs increases with frequency with a slope result- 
ing from a balance between the excitation and damping 
rates. Nevertheless, modes of angular degree i = 1 exhibit 
a singular behavior, i.e. a maximum at v ~ 60 /xHz. This is 
due to the variation of the slope in the excitation rates (see 
Fig. 5). In terms of amplitudes, the maximum is found to 
be w 5mm/s for ^ = 1 at w 60 /iHz, which corresponds 
to the mode with radial order |n| = 10. It is important 
to stress that the velocities shown in Fig. 10, are intrin- 
sic values of the modulus that must not be confused with 
the apparent surface velocities (see Sect. 5), which are the 
values that can be compared with observed ones. 



4.2. Comparison with previous estimations 

The theoretical intrinsic velocities obtained in the present 
work must be compared to previous estimations based on 
the same assumption that modes are stochastically excited 
by turbulent convection. All works cited in the next sec- 
tions deal with intrinsic velocities, i.e. ones not corrected 
for visibility effects. 



4.2.1. Estimation based on the equipartition of energy 

The first estimation of g-mode amplitudes was performed 
by Gough (1985), who found a maximum of velocity of 
about 0.5mms~^ for the i = 1 mode at « 100 /xHz. 
Gough (1985) used the principle of equipartition of energy, 
which consists in equating the mode energy {£) with the 
kinetic energy of resonant eddies whose lifetimes are close 
to the modal period. This "principle" has been theoreti- 
cally justified forp modes, by Goldreich & Keeley (1977b) 
assuming that the modes are damped by eddy viscosity. 
They found that the modal energy to be inversely propor- 
tional to the damping rate, 77, and proportional to an in- 
tegral involving the term E\ vx X where E\ = (1/2) m\ 
is the kinetic energy of an eddy with size A, velocity v\, 
and mass m\ = pX^ (see Eq. (46) of Goldreich & Keeley 
1977b). Using a solar model, they show that the damping 
rates of solar p modes are dominated by turbulent viscos- 
ity and that the damping rates are accordingly propor- 
tional to the eddy-viscosity, that is, rj (x vxX (see Eq. (6) 
of Goldreich & Keeley 1977a). Hence, after some simplify- 
ing manipulations, Goldreich & Keeley (1977b) found the 



10 



Belkacem et al.: Stochastic excitation of nonradial modes II. Are solar asymtotic gravity modes detectable? 



modal energy to be (see their Eq. (52)) 
£ « 0.26 = 0.13 mA v\ 



(23) 



This principle then was used by Christcnscn-Dalsgaard 
& Frandsen (1983) for p modes and Gough (1985) for 
solar g modes. However, as mentioned above, the result 
strongly depends on the way the modes arc damped, and 
for asymptotic g modes there is no evidence that this ap- 
proach can be used and in particular, as shown in this 
work, if the damping is dominated by radiative losses. 

4.2.2. Kumar et al. (1996)'s formalism 

Another study was performed by Kumar et al. (1996), 
which was motivated by a claim of g-mode detection in 
the solar wind (Thomson et al. 1995). Computations were 
performed using the Goldreich et al. (1994) formalism; 
both turbulent and radiative contributions to the damp- 
ing rates were included as derived by (Goldreich & Kumar 
1991) who obtained mode lifetimes around 10^ yrs. This 
is not so far from our results (see Fig 9). Kumar et al. 
(1996) found that the theoretical (i.e. not corrected for 
visibility factors) surface velocity is around 10~^cms~"'^ 
near v — 200 /iHz for I ~ \ modes. However, as shown 
in Sect. 3, the results for this frequency range are very 
sensitive to the convective flux perturbation in the damp- 
ing rate calculations. Thus, we do not discuss the result 
obtained for those frequencies. 

More interesting for our study, Kumar et al. (1996) also 
found very low velocities (10~^ mms"-'^) for v < 100 /^Hz. 
This is significantly lower than what we find. However, 
the efficiency of the excitation strongly depends on how 
the eddies and the waves are temporally-correlated. As 
already explained in Sect. 2.1, the way the eddy-time cor- 
relation function is modeled is crucial since it leads to 
major differences between, for instance, a Gaussian and 
a Lorentzian modeling. The Goldreich & Keeley (1977b) 
approach, from which Kumar et al. (1996)'s formulation is 
derived, implicitly assumes that the time-correlation be- 
tween the eddies is Gaussian. The present work (as ex- 
plained in Sect. 2.1) assumes a Lorentzian for the time 
correlation function Xk, which results in v ~ 3 mm/s in 
amplitude ior £ = 1 mode at ;y w 60 fiHz (Sect. 4.1). 

We performed the same computation but now assum- 
ing Xk to be Gaussian (Eq. (7)) and using a Kolmogorov 
spectrum as in Kumar et al. (1996). In that case (see 
Fig. 10), we find velocities of the order of 10~^TOms~^ 
for £ = 1 which agree with the result of Kumar et al. 
(1996), which is significantly lower than when assuming a 
Lorentzian. 



5. Apparent surface velocities 

We denote as disk-integrated apparent velocities the values 
of amplitudes that take both geometrical and limb dark- 
ening effects into account. Contrary to solar p modes, one 
cannot neglect the horizontal component of ^ compared 
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Fig. 10. Top: Theoretical intrinsic surface velocities of g- 
modes of degree ^ = 1,2,3 as a function of the oscillation 
frequency in /xHz, computed as described in Sect. 4.1 using a 
Lorentzian Xk- Bottom: Surface velocities of gravity modes of 
angular degree 1=1 and £ = 2 computed using a Gaussian Xk 
and a Kolmogorov spectrum to reproduce the results of Kumar 
et al. (1996). 







1=1 0.117 
1=2 0.346 
1=?, 0.06 



0.675 
0.107 
0.164 



0.437 
0.0552 



0.184 



Table 1. Values of the visibility coefficient a™ of the radial 
component of the velocity, corresponding to an inclination an- 
gle of Oo = 83°. 



to the vertical one. The observed velocity [Vobs] is given 
by the apparent surface velocity < |T4pp(r, t)p >^/^(see 
Appendix. C) evaluated at the observed line formation 
height h: 



Vr. 



■bs 



p 



1/2 



{a'T Vrih) + (3T vnih)) (24) 



where a™ and /3™ are the visibility factors defined in 
Appendix. C. 
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m 1 2 3 

i^l 0.094 0.540 

£=2 0.833 0.258 1.053 

1=3 0.291 0.649 0.268 0.892 

Table 2. Values of the visibility coefBcient /3™ of the radial 
component of the velocity, corresponding to an inclination an- 
gle of Bo = 83°. 
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Fig. 11. Apparent surface velocities for g modes of degree 
£ = 1,2,3 as a function of the oscillation frequency in /iHz 
(visibility factors are taken into account). 



In Appendix C we follow the procedure first derived 
by Dziembowski (1977) and for asymptotic g modes by 
Berthomieu & Provost (1990). We use a quadratic limb- 
darkening law following Ulrich et al. (2000) for the Sun 
with an angle between the rotation axis and the Equator 
of 83 degrees. As mentioned above, the apparent velocities 
are evaluated at the level h, i.e. the height above the pho- 
tosphere where oscillations arc measured. Then h is set 
so as to correspond to the SoHO/GOLF measurements 
that use the NaDl and D2 spectral lines, formed at the 
optical depth r = 5.10"'' (see Bruls & Rutten 1992). The 
results are presented in Tables 1 and 2 for angular degrees 
i = 1,2,3. 

Figure 11 displays the apparent velocities for modes 
£ = 1,2,3 and £ = m. For a given angular degree, the 
azimuthal order degree is chosen such that the apparent 
velocity is maximal. The velocities of the m = modes are 
strongly attenuated by the visibility effects, while the m = 
£ modes are less sensitive to them. For £ = 1 modes, the 
amplitudes are divided by a factor of two with respect to 
the intrinsic velocities, while the £ = 2,3 mode velocities 
remain roughly the same. Consequently, our calculations 
show that both the £ = 1 and £ = 2 {m = £) are the 
most probable candidates for detection with amplitudes 
w 3 mm s~^. 



5.1. Detectability of g modes; only a matter of time 

To compare our calculated apparent velocities with ob- 
servations, we used data from the GOLF spectrome- 
ter (Gabriel et al. 2002) onboard the SOHO platform, 
which performed Doppler-like measurements on the disk- 
integrated velocity of the Sun, using the Na D lines. We 
used here a series of 3080 days to estimate the background 
noise level and compare it to the apparent velocities de- 
termined in this work. 

A first possible approach is to use some analytical 
and statistical calculations such as the ones developed by 
Appourchaux et al. (2000) (Eq. 10). Once a length of ob- 
servation T (in units of lO^s), a frequency range Au (in 
//Hz), and a level of confidence pdet are set, this gives the 
corresponding signal-to-noisc ratio 

^ ln(r) + ln(Ai.) - ln(l - p,et) , (25) 

< s > 

where s^i^t is the power of the signal to be detected, and 
< s > the local power of the noise. This means that any 
peak in the frequency range Aiy above this ratio has a 
probability pdet of not being due to noise. Choosing a fre- 
quency range of Av = 10/xHz centered on the frequency 
of the highest expected velocities (60/iHz) sets the back- 
ground level at w 500 (ms~^)^/Hz. Equation (25) gives an 
amplitude of 5.2 mm for a detection with a confidence 
level Pdet of 90% for 15 years of observation, 4.6 mms^^ 
for 20 years, and 3.8 mms~^ for 30 years. 

However, this approach has to be repeated for each 
mode (with its own proper noise level) to have a global 
view of detection possibilities. To do so, we used sim- 
ulations. Again relying on the GOLF data to estimate 
the noise spectrum, we simulated synthetic data including 
noise and g modes with the apparent velocities as above 
(and with random phases). Several durations of observa- 
tion were simulated, from 10 to 30 years. A hundred sim- 
ulations were performed in each case. The noise level is 
estimated locally and so is the detection level, following 
Eq. (25), on the frequency range [30/fHz,100/iHz]. Thus, 
with a confidence level of 90% and with 7 independent 
subsets of lOfiKz, noise is expected to show no peak above 
the global detection level with a probabihty of 48%, and to 
show 1 peak above the global detection level with a proba- 
bility of 32% (and even 2 peaks in 12% of the realizations). 
Table 3 lists the average (over 100 simulations) number 
of peaks detected above the detection level for different 
observation durations. These simulations were performed 
using amplitudes A„iax assuming three different cases: 

• Case 1: we assumed for Amax the apparent surface 
velocity amplitudes calculated above, A. 

Due to uncertainties in the theoretical modeling (as 
discussed in Sect.6), we also assume: 

• Case 2: that amplitudes are larger than the ampli- 
tudes estimated above by 50% i.e. A,nax = 1.5 A 

• Case 3: that amplitudes are larger than the ampli- 
tudes estimated above by a factor 2 i.e. A„iax ^ 2 A 

Cases 1 and 3 are the two limits of this exercise. The 
number of detected peaks in Case 3 shows that the pre- 
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Fig. 12. Simulated spectrum for an observation length of 30 
years, in Case 1. The dashed line indicates the level of de- 
tection (see text and Eq. 25). The vertical lines indicate the 
frequencies of the simulated modes. Here, only one mode is 
above the detection level. Bottom: Average number and stan- 
dard deviation (from the 100 realizations) of modes detected 
in simulations versus the length of observation for the three 
cases (see text). Above the hashed region (less than four peaks 
detected), one can consider the detection to be unambiguous. 
The upper curve corresponds to the case 3 (Amax = 2^4), the 
middle one to case 2 (ylmax = 1.5A), and the lower one to case 

1 (^max = A) • 

dieted amplitudes cannot be overestimated by a factor of 
two, because in this case, the solar g modes would have 
already been detected without doubt. Case 1 sets a lower 
limit, because in this case, even with longer (30 years) ob- 
servation, g modes would not be detected. Case 2 shows 
that if real solar amplitudes are just a few tens of percent 
higher than the present estimations, then g modes could 
be detected no doubt after say 15 to 20 years of observa- 
tion (to be compared to the present status of observation: 
12 years). The results are summarized in Figs. 12. 

We must stress that, apart from visibility effects and 
height of line formation, we took no other instrumental 
effects on the apparent amplitude determination into ac- 
count, because they depend on the instrument. The im- 
pact is probably a decrease in the measured amplitudes 
compared to the apparent amplitudes as computed here. 



Amax 10 years 15 years 20 years 25 years 30 years 

' A 08 Te OS O O 

1.5A 1.4 2.9 4.5 6.5 8.8 

2A 4.6 8.5 13.4 20.0 21.~ 

Table 3. Number of peaks above the detection level in the 
simulated power spectra versus the duration of observation in 
three cases. In the simulated signal, the modes are given an 
amplitude Amax- The 3 cases respectively correspond to Amax 
being the apparent amplitudes A readily stemmed from our 
calculation, Amax = 1.5^ and Amax = 2j4. The last two cases 
take into account that uncertainties in the modeling globally 
tend to underestimate the amplitudes as discussed in Sect. 6. 



This does not change the above conclusion for Case 1. We 
expect that the instrumental uncertainty is less than the 
theoretical uncertainties discussed in Sect. 6 below, which 
led to Case 2 and 3. 

6. Discussion 

In Sect. 5 above, we explained why estimates of g- 
mode amplitudes obtained by previous authors differ from 
each other by orders of magnitude (Christensen-Dalsgaard 
2002b). We propose an improved modeling based on the 
input of 3D numerical simulations and on a formalism that 
had successfully reproduced the observations for p modes 
(Belkacem et al. 2006). Nevertheless, several approxima- 
tions remain, and they lead to uncertainties that can reach 
a factor two in the estimation of g-mode apparent veloci- 
ties (overestimation). We next discuss the most important 
ones. 

6.1. Equilibrium model: description of convection 

Convection is implemented in our equilibrium models 
according to the classical Bohm-Vitense mixing-length 
(MLT) formalism (see Samadi et al. 2006, for details). 

6.1.1. Convective velocities 

Values of the MLT convective velocity, u, are by far the 
most important contributions to mode amplitude uncer- 
tainties because the mode surface velocities depend on u'^ 
. First, we verified that a non-local description of turbu- 
lence does not modify the convective velocities by more 
than a few per cent except near the uppermost part of 
the convection zone, which does not play any role here. 
Second, we compared the rms velocities from the 3D nu- 
merical simulation with MLT velocities to estimate of the 
uncertainties. The MLT underestimates the velocity, rel- 
ative to the more realistic numerical simulation (far from 
the boundaries). Indeed, it comes from the negative kinetic 
energy flux that results in a larger enthalpy flux in order 
to carry the solar flux to the surface. A direct consequence 
is that in 3-D simulations the velocities are higher than 
the ones computed by MLT by a factor of about 50%. 
This may in turn result in a possible underestimation of 
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the amplitudes of the modes by a factor 2, when, as here, 
MLT is used to estimate the velocities. 

6.1.2. Anisotropy 

The value for the velocity anisotropy, which is the ratio 
between the square of the vertical velocity to the square of 
the rms velocity parameter, $, is derived from the MLT: 
its value is 2. However, this is not fully consistent since we 
assume, in the excitation model, isotropic turbulence (i.e. 
4> = 3). Nevertheless, increasing the value of $ from two to 
three results in an increase of only 15 % in the mode sur- 
face velocities. This is lower than the uncertainties coming 
from Xk (see Sect. 6.2). 

6.1.3. Turbulent pressure 

Our solar equilibrium model does not include turbulent 
pressure. However, unlike p modes, low-frequency (high 
radial order) gravity modes, i.e. those considered in this 
work, are only slightly affected by turbulent pressure. The 
reason is that such modes are excited in the deepest lay- 
ers of the convection zone, i.e. between r = 0.7 Rq and 
r = 0.9 Rq where turbulent pressure has little influence 
on the equilibrium structure since the ratio of the turbu- 
lent pressure to the gas pressure increases with the radius. 

6.2. Stochastic excitation: the role of the eddy-time 
correlation function 

A Gaussian function is commonly used to describe the 
frequency dependence of the turbulent kinetic energy 
spectrum, Xfc j (e-g-, Samadi & Goupil 2001; Chaplin 
et al. 2005). However, Samadi et al. (2003a) show that, 
for p modes, a Lorentzian function represents the re- 
sults obtained using 3D numerical simulations better. 
Furthermore, the latter function yields a theoretical mod- 
eling in accordance with observations, while using a 
Gaussian function fails (Samadi et al. 2003b). This led 
us to investigate Xki^^) for g modes. We find that differ- 
ent choices of the functional form for Xkii-^) result in order 
of magnitude differences for the mode amplitudes. 

Uncertainties inherent in the eddy-time correlation 
function are related to the value of the A parameter 
(Sect. 2.2.1) and to the contribution of low frequency com- 
ponents in the 3D simulation. As a rough estimate, de- 
creasing A from 3 to 2 leads to an increase of 20 % for 
the surface velocity. Figure 3 shows that low-frequency 
components in the turbulent kinetic energy spectrum 
are better-fitted using a Gaussian function. However, the 
source of such low-frequency components remains unclear, 
because they can originate from rotation; in particular, it 
is not clear whether they must be taken into account when 
estimating the mode excitation rates. By removing those 
contributions, the resulting surface velocities decrease by 
around 25 %. 



6.3. Mode damping: the convection-pulsation coupling 

Last but not least, modeling damping rates of damped, 
stochastically excited modes remains one of the most chal- 
lenging issues. The strong coupling between convection 
and oscillation in solar-like stars makes the problem diffi- 
cult and still unsolved, since all approaches developed so 
far failed to reproduce the solar damping rates without 
the use of unconstrained free parameters (e.g., Dupret et 
al. 2005, Houdek 2006). Such descriptions fail to correctly 
describe the interaction between convection and oscilla- 
tions when both are strongly coupled, i.e. when the char- 
acteristic times associated with the convective motions are 
the same order of magnitude as the oscillation periods. 
This explains why we do not use an extrapolation based 
on a fit of p mode damping rates, but instead consider 
a frequency domain in which the damping is dominated 
by radiative contributions. A reliable computation of the 
damping rates at higher frequencies, beyond this paper's 
scope, would require a sophisticated analytical or semi- 
analytical theory of the convection-oscillation interaction, 
which will not be limited to the first order in the convec- 
tive fluctuations and which will take the contribution of 
different spatial scales into account. 

7. Conclusions 

We performed a theoretical computation of the surface 
oscillation velocities of asymptotic gravity modes. This 
calculation requires knowing excitation rates, which were 
obtained as described in Belkacem et al. (2008) with in- 
put from 3D numerical simulations of the solar convec- 
tive zone (Miesch et al. 2008). Damping rates, t] , are 
also needed. As mentioned in Sect. 6, we restricted our 
investigation to the frequency domain for which r] is dom- 
inated by radiative contributions {i.e. v S [20; 110]/xHz). 
For higher frequencies, the coupling between convection 
and oscillation becomes dominant, making the theoreti- 
cal predictions doubtful. For asymptotic f/-modes, we find 
that damping rates are dominated by the modulation of 
the radial component of the radiative flux by the oscilla- 
tion. In particular for the I = \ mode near « 60 /xHz, r\ 
is around 10^^/iHz, then the mode life time is « 3.10^yrs. 
Maximum velocity amplitude at the photosphere arises for 
this same mode and is found at the level of 3 mms~^ (see 
Fig. 11). Modes with higher values of the angular degree 
I present smaller amplitudes since the damping is propor- 
tional to . 

Amplitudes found in the present work are orders of 
magnitude larger than those from previous works, which 
themselves showed a large dispersion between their respec- 
tive results. In one of these previous works, the estimation 
was based on an equipartition principle derived from the 
work of Goldreich & Keeley (1977a,b) and designed for 
p modes. Its use for asymptotic g modes is not adapted 
as the damping rates of these modes are not dominated 
by turbulent viscosity. Kumar et al. (1996) have carried 
another investigation of 5 mode amplitudes, and its calcu- 
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lation is rather close to our modeling. Most of the quanti- 
tative disagreement with our result lies in the use of a dif- 
ferent eddy-time correlation function. Kumar et al. (1996) 
assumed a Gaussian function as is commonly used. Our 
choice relies on results from 3D simulations and is closer 
to a Lorenztian function. 

Taking visibility factors, as well as the limb-darkening, 
into account we finally found that the maximum of ap- 
parent surface velocities of asymptotic ij-modes is « 3 
mms~^ for £ = 1 at « 60/iHz an £ = 2 at « 100 /xHz. 
Due to uncertainties in the theoretical modeling, ampli- 
tudes at maximum, i.e. for ^ = 1 at 60 //Hz, can range 
from 3 to 6 mms^^. By performing numerical simulations 
of power spectra, it is shown that, with amplitudes of 6 
mms~^, the modes would have been already detected by 
the GOLF instrument, while in the case of an amplitude of 
3 mms~^ the g modes would remain undetected even with 
30 years of observations. The theoretical amplitudes found 
in this work are then close to the actual observational 
limit. When detected, the amplitude detection threshold 
of these modes will, for instance, establish a strict upper 
limit to the convective velocities in the Sun. 
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Appendix A: Energy equation near the center 

For the full non-adiabatic computation of i^-mode damp- 
ing rates, much care must be given to the solution of the 
energy equation near the center of the Sun for the modes 
of angular degree £ = 1. We give in Eqs. (A.l) and (A. 2) 
the perturbed energy and transfer equations in a purely 
radiative zone: 
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(A.2) 



dT/dr dr 

The radial (first term of Eq. (A.l)) and transverse 
parts (last term of Eq. (A.l)) of the perturbed flux diver- 
gence arc both singular at the center. But this singularity 
is lifted when the two terms are joined and an appropriate 
change of variables is carried out: 
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All of these variables and quantities are regular at the 
solar center, where the perturbed energy equation takes 
the form 
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(A.4) 



For a precise solution of the non-adiabatic problem by 
a finite difference method, it is crucial to use a discrete 
scheme that tends continuously towards Eq. (A.4) at the 
center. If not, the eigenfunctions diverge towards the cen- 
ter; in the particular case of the solar g modes, this can 
lead to an overestimate of the damping rates by a factor 
of about 2. 
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Appendix B: Computation of the kinetic energy 
spectrum from the ASH code 

The ASH code solves the hydrodynamical equations in 
spherical coordinates (r, 6, 0). Each component of the ve- 
locity is decomposed in terms of spherical harmonics as 



(B.l) 



where p ~ 
defined as 



The spherical harmonic Y;_^ 



y/"(0,0)=iV,,„Fr(cos0) e^™* 



(B.2) 



where is the associated Legendre function, and the 
normalization constant iVj.m 



is chosen such that 
dnYi,m{0, 
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(B.4) 



where dil = sm OdOdcf). 

The kinetic energy spectrum that is averaged over 
time and the solid angle is defined following Samadi et al. 
(2003b) as 



{yi,7n,p)y 



(B.5) 



where ((.)) refers to time average. As in Samadi et al. 
(2003b), density does not enter into the definition of the 
kinetic energy spectrum. Indeed, the Samadi & Goupil 
(2001)' formalism assumes a homogeneous turbulence. 
This assumption is justified when the turbulent Mach 
number is low. This is the case in most parts of the con- 
vective zone except at the top of convective region. 

The mean kinetic energy spectrum, E[l, r), verifies the 
relation 

J2E{l,r)^lu\r) (B.6) 
I 

where u{r) is the root mean square velocity at the radius 
r. 

Following Samadi et al. (2003a), we also define a ki- 
netic energy spectrum as a function of frequency (u) and 
averaged over the solid angle, E{1, v^r) such that 



\Vp{v,r,e,4>) (B.7) 



where Vp{v^r^9,(j}) is the time Fourier transform of 
Vp{t,r,e,(j)) - {Vp). Using Eqs.(B.l) and (B.4), Eq. (B.7) 
yields: 



where V;,rn,p('^j^) is the time Fourier transform of 
yi,m,p{i,'r) — {Vi^„L,p)- As in Samadi et al. (2003a), we de- 
compose E{1, v, r) as 



E{l,y,r)^E{lr) xM,r) 



(B.9) 



where the fimction xi{^j satisfies the normalization con- 
dition 



+ 00 



dvxi{y,r) = l. (B.IO) 
According to the Parseval-Plancherel relation, one has 

/ + 00 1 
du Eil, v,r) = Y, E{1. r)^- u'ir) , (B.ll) 
I ^ 

We consider a short time series of duration « 4.68 days 
with a sampling time of 800 seconds. Accordingly the 
Nyquist frequency is w 1 mHz and the frequency resolu- 
tion reachs « 2.5 /iHz. In addition, we use a longtime se- 
ries of duration « 45.83 days with a sampling time of 4 10** 
seconds that permits us to get Xk at very low frequencies. 
In practice, E(X) is derived from Eq. (B.5) and is directly 
implemented into Eq. (1), while Xfc(j^) inferred from the 
simulation is computed using Eqs. (B.9) and (B.5). 

By using E(X) from the numerical simulation, we as- 
sume a planparallel approximation {E(k)dk = E{l)dl) 
since the maximum of the kinetic energy spectrum occurs 
on scales ranging between / w 20 and Z w 40. 

Appendix C: Visibility factors 

Visibility factors have been computed first by 
Dziembowski (1977). Berthomieu & Provost (1990) 
studied the case of g modes which, for convenience, we 
recall below in our own notation. We denote the spherical 
coordinate system in the observer's frame by (r, 0, (p) 
where r = corresponds to the center of the star and the 
9 = axis coincides with the observer's direction. At a 
surface point (r, 0, 0), the unit vector directed toward the 
observer is n = cos^e^ — smOeg. The apparent surface 
velocity is obtained as 



^app (^: ^) 



Jh{^i) {v{r,t)-n)dn 
/ h{fi) dn 



(C.l) 



where v{r,t) is the intrinsic mode velocity and h{fi) the 
limb-darkening function, which is normalized such that: 



fi h{^)dii = 1 . 



(C.2) 



m,p 



To first order in linearized quantities in Eq. (C.l), the 
effect of the distorted surface is neglected, and dfl = 

sin OdOdcj) is the solid angle around the direction of the 
observer n with R the stellar radius. 

For slow rotation, the oscillation velocity can be de- 
scribed in a pulsation frame with a single spherical har- 
monic. The coordinate system (r, Q, $) in the pulsation 
frame is chosen such that the pulsation polar axis coin- 
cides with the rotation polar axis. The velocity vector at 
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a level r in the atmosphere of the star for a mode with As emphasized by Dziembowski (1977), only the qf^ q 

given i,m and pulsation frequency loq can then be writ- coefficients survive the 4> integration in Eq. (C.l). Its ex- 

tcn with no loss of generality as pression is 

v{r, t) i a{t) ^0 |(r) e^""' + c.c. (C.3) 9^,0(00, $0) = Pr(cos 60) e^™*° . (C.ll) 

where c.c. means complex conjugate and with the displace- The angle 60 between the observer and the rotation axis 

ment eigenvector defined as often denoted i. Integration over the solid angle leads 



to 



iir) = &{r) Yne, $) Br + VHYne, $) (C.4) 

with 



v„^(o±J_±) (C5) Ur{r) ^i''h{^^)p,{^i)d^i + 

" ^ '99' sine 9$^' ^ ' ^ Jo 

The dimensionless complex velocity amplitude a«(t) ^h{r) ( H h{ii) (1 — ^^) ^^^^d/i^ (C12) 
is assumed to be a slowly varying function of time for M 

a damped stochastically excited mode (Samadi & Goupil Finally, using properties of spherical harmonics, one ob- 
2001; Samadi et al. 2003b; Belkacem et al. 2008). The ^^^^^ 
theoretical expression is given by 

pi{n){i{r) ■ n)dVL ^ \ (c ( \ ( \ \ 

< laWP >= A (C.6) IH,)dn ^^'^^'^ + 

Vl^o (C.13) 
where the power P is defined in Eq. (1), / is the mode "^^^^^ ^^"'^ defined 
inertia, 7] the damping rate and <> represents a statistical ri 

average, or cquivalcntly here a time average. = / d/x/i h{fi) Pi{ii) (C.14) 

To obtain the apparent velocity from Eq. (C.l) using " ^ 

Eqs. (C.3) and (C.4), one must compute the scalar prod- ^Ji = £ [ d/x ^ h{fi) {Pi_i - ^Pg) (C.15) 
uct: $,{r) ■ n. Jo 

|(r) . n = Ur) ¥^0, <J>) (e. • n) + ^(r) (V^l^r ' ") • w . 

(C.7) Mm)= ■ (C.16) 

Jo h[fi)dfi 

A change in coordinate system shows that • n = cos 6* 

and Collecting Eq. (C.3) and Eq. (C.13), the apparent ve- 

VhYP (6, $) . n = - sin e^^HMl lo^ity is then given by 

06 ^ 

We use the spherical harmonics as defined in Eq. (B.2) and Vapp{r, t) ~ - a{t) ujq Ni^m Pi(co& 60) (C.17) 



the following property 

I 

Pr(cose)e™*= 4>n'(0o,$o)P^'(cos0)e™'^ 

m' = -t 

(C.8) 

which for convenience, we use under the form 



2 

(CIS) 

We assume a quadratic limb-darkening law of the form 

h{^)^l + ci +C2 (C.19) 

^ where X = I — Ci={i_2.3} are the associated limb- 

y™(Q $) — TVf (60 $0) ^""'(cos^) e"" '^ darkening coefficients, which respective values are —0.466, 

^ ' '™ "T"^ „ ™'™ ' ^ —0.06 and —0.29 for the NaDl spectral line, as derived 

(C 9) ^■^ Ulrich et al. (2000). We find that our conclusion de- 
pends neither on the adopted limb-darkening law nor on 

where Ni^,n is defined in Eq. (B.3), and (60, $0) are the the Hmb-darkening coefficients, results in accordance with 

coordinates of the line-of-sight direction in the pulsation Berthomieu & Provost (1990). 

frame. The scalar product Eq. (C.7) becomes Using Eq. (C.6), the rms velocity is obtained as: 



$(r) • n - m,ra 'Z™,™'(eo,$o) e™> (< |T4pp(r,t)P >y/^ = ( m,m \Piicoseo)\ 

(C.20) 
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which we finally rewrite as 

«IWM)P»V^=(^)''' (C.21) 
X \vr{r) aT + Vh{r) P^l, 

where we have defined 

aT = Ne^m PKcosGo) ui (C.22) 
f3T = Ne,„^ Pi{coseo)wi, (C.23) 

and 

Vr{r) = ujQ ^rir) ; Vh{r) = ojq ^hir) ■ (C.24) 



